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ABSTRACT 


This  paper  covers  the  design  of  a  dynamic  positioning  system  for  the  U.S.  Coast  Guard 
WLB  "IRIS"  class  buoy  tender.  The  control  system  design  is  based  upon  optimal  con¬ 
trol  theory  with  estimates  of  position  and  heading  provided  by  a  steady  state  Kalman 
filter.  Sea  current  estimates  are  provided  by  a  Kalman  filter  predictor  based  upon  the 
innovations  process.  The  vessel  and  dynamic  positioning  system  are  simulated  using  the 
Dynamic  Simulation  Language  (DSL). 
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THESIS  DISCLAIMER 


The  reader  is  cautioned  that  computer  programs  developed  in  this  research  may  not 
have  been  exercised  lor  all  cases  of  interest  While  every  effort  has  been  made,  within 
the  time  available,  to  ensure  that  the  programs  are  free  of  computational  and  logic  er¬ 
rors,  they  cannot  be  considered  validated.  Any  application  of  these  programs  without 
additional  verification  is  at  the  risk  of  the  user. 
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I.  INTRODUCTION  TO  DYNAMIC  POSITIONING  SYSTEMS 

A.  DEFINITION 

A  dynamic  positioning  system  (DPS)  is  an  automatic  control  system  that  will  drive 
a  vessel  to  a  selected  position  and  maintain  that  position  and  selected  heading  using 
computer  control  of  propulsive  thrusters.  A  DPS  includes  one  or  more  position  and 
heading  measurement  systems  and  a  computer-directed  control  system.  The  system 
should  keep  the  vessel  within  its  specified  position  and  heading  constraints  with  optimal 
use  of  propulsors,  reducing  the  wear  on  that  machinery. 

B.  ADVANTAGES  AND  DISADVANTAGES 

The  use  of  DPS  is  not  new.  The  first  dynamic  positioning  systems  appeared  in  the 
late  1960's  on  coring,  cable  laying  and  surface-support  ships  for  underwater  work.  The 
success  of  these  initial  DP  vessels  spread  the  application  of  the  systems  into  offshore  oil 
industry,  dredging,  precision  dumping,  pipe  laying,  and  floating  hotels. 

There  are  several  advantages  of  dynamic  positioning  over  anchor  deployment.  It  is 
usable  in  all  water  depths  and  not  constrained  by  the  length  of  its  anchor  chain.  In 
contrast  to  anchoring,  there  is  virtually  no  set-up  time  required  to  begin  stationkccping. 
In  order  to  reach  and  attain  a  precision  position  by  anchoring,  time  and  personnel  in¬ 
tensive  anchoring  details  must  be  set.  A  DPS  allows  a  vessel  to  attain  virtually  any  po¬ 
sition  and  heading  automatically.  Due  to  its  dynamic  nature,  a  DPS  vessel  can  work  in 
close  quarters  with  fixed  platforms  or  anchored  vessels.  A  dynamic  positioning  system 
is  designed  to  operate  in  adverse  weather  conditions 

The  disadvantages  of  DPS  include  high  initial  capital  investment,  high  fuel  costs 
(over  anchoring),  high  maintenance  due  to  the  continuous  activation  of  thrusters,  and 
increased  manpower  requirements  for  maintenance. 

Dynamic  positioning  is  an  effective  means  of  stationkeeping  when  water  depth  is  too 
great  for  anchoring,  on-station  time  will  be  brief,  the  vessel  is  required  to  be  on  station 
in  all  weather  conditions,  and  accuracy  in  positioning  is  more  important  than  cost. 

C.  EQUIPMENT  FOR  DYNAMIC  POSITIONING  SYSTEMS 

The  equipment  needed  for  a  dynamic  positioning  system  includes  a  wind  sensor 
(velocity  and  direction),  gyrocompass,  and  a  position  reference  sensor. 
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There  are  several  position  seniors  that  are  employed  in  DPS  designs.  They  include 
sensors  that  determine  position,  velocity  and/or  acceleration  in  some  common  reference 
system.  The  basic  position  sensor  types  are  taut  wire,  optical,  acoustic  beacon,  and  ra¬ 
dio  svstems. 

'I'he  taut  wire  scheme  works  by  lowering  a  weight  to  the  sea  floor  by  wire  rope.  The 
wire  is  held  in  constant  tension,  usually  by  a  constant-tension  winch,  and  any  angle  de¬ 
viations  from  vertical  are  corrected  by  thrusters. 

Optical  positioning  systems  are  of  either  active  or  passive  design  and  make  use  of 
laser  technology.  Triangulation  of  sextant  angles  provides  position  information.  How¬ 
ever,  sextants  are  manually  trained  by  personnel. 

Acoustic  beacon  is  the  most  common  positioning  system  in  the  offshore  industry. 
This  system  requires  a  sea  floor  beacon  and  hull-mounted  subsurface  hydrophones. 

In  general,  standard  radionavigation  systems  such  as  Loran-C,  OMEGA,  and  Tansit 
SATNAV  do  not  provide  sufficient  accuracy  for  DP  systems.  The  approximate  accura¬ 
cies  of  these  systems  are: 

•  Loran-C  -  460  meters, 

•  OMEGA  -  500  meters, 

•  SATNAV  -  25  meters. 

[Ref.  lr:p.  8-18) 

Accuracies  of  those  magnitudes  are  completely  unacceptable  for  DP  systems.  The 
NAVSTAR  Differential  GPS,  however,  is  able  to  provide  accurate  position  information 
to  within  3  meters  anywhere  on  the  globe.  In  addition,  this  system  allows  vessels  to 
measure  velocity  using  Doppler  to  accuracies  of  about  0.1  meters  per  second.  [Ref.  1  : 

p.  8-18). 

For  a  ship  to  be  dynamically  positioned  it  must  use  thrusters,  sensors,  computer, 
and  a  control  and  display  system.  The  thrusters  must  be  able  to  provide  fore  and  aft, 
atlnvartships,  and  moment  control.  Sensors  provide  physical  information  including  po¬ 
sition,  heading,  and  wind  conditions.  The  computer  uses  sensor  information  and  calcu¬ 
lates  required  thruster  commands  that  will  position  the  vessel  at  th;  reference  position 
and  heading.  Controls  and  displays  are  required  to  allow  the  operator  to  monitor  and 
control  the  DP  system. 

The  purpose  of  this  thesis  is  to  design  and  simulate  a  dynamic  positioning  system 
for  a  li.S.  Coast  Guard  buoy  tender.  The  dynamic  positioning  system  is  designed  to  be 
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fully  automated,  performing  its  functions  in  varying  sea  states,  wind,  and  current  con¬ 
ditions  with  its  only  operator  input  being  the  referenced  position. 
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II.  BACKGROUND 


The  U.S.  Coast  Guard  has  begun  the  acquisition  process  for  seagoing  and  coastal  « 

buoy  tenders  to  replace  the  aging  157  foot  WLM  class  and  the  180  foot  WLB  class 
tenders.  Sec  Figure  1  for  a  comparison  of  the  old  WLB  fleet  with  its  successor  [Ref.  2: 
p.7).  One  of  the  requirements  for  the  replacement  vessels  is  a  dynamic  positioning  sys¬ 
tem.  The  dynamically  positioned  vessel  must  be  able  to  approach,  maneuver,  and 
maintain  poruicn,  within  a  circle  of  a  10  meter  radius  over  a  fixed  point  on  the  earth, 
alongside  floating  aids  to  navigation  in  restricted  channels  of  shallow,  18  foot  deep  bays, 
estuaries,  and  rivers.  The  system  is  required  to  maintain  position  in  a  30  knot  wind  of 
any  aspect  and  5  knot  current  within  10  degrees  of  ship's  head  [Ref.  3  :  p.  15  ].  This 
DPS  design  endeavors  to  take  a  modified  existing  U.S.  Coast  Guard  (USCG)  buoy 
tender  using  a  minimum  of  additional  control  and  measurement  inputs  in  order  to  reduce 
the  potential  implementation  costs.  This  initial  design  is  a  180  foot  WLB  "IRIS"  class 
buoy  tender.  The  only  modifications  to  existing  vessels  are  the  installation  of  a  stern 
thruster  and  Global  Positioning  System  (GPS).  The  GPS  provides  real-time  positioning 
information  accurate  to  within  three  meters  anywhere  on  the  globe  [Ref.  1  :  p.  8-18  ). 

The  GPS  is  a  requirement  of  the  replacement  buoy  tender  fleet.  The  available  controls 
for  this  design  arc  200  S11P  bow  and  stern  thrusters  and  1170  SUP  single  screw  propul¬ 
sion.  System  sensor  measurements  are  GPS,  gyrocompass,  and  wind  inputs.  The  mul¬ 
tiple  mission  requirements  of  the  buoy  tender,  as  with  all  Coast  Guard  cutters,  placed 
some  constraints  on  the  DPS  design  from  the  outset.  The  buoy  tender  is  occasionally 
called  upon  to  perform  some  harbor  icebreaking  duties.  It  is  considered  an  ice-capable 
vessel,  able  to  break  up  to  one  foot  of  ice  at  three  knots  [Ref.  3  :  p.  15  ).  Ice  operations 
preclude  the  use  of  subsurface  position  sensing  devices  such  as  hydrophones  and  taut 
wire  schemes.  In  order  to  clarify  the  reasoning  behind  some  of  the  system  constraints, 
a  brief  qualitative  description  of  buoy  tending  operations  follows. 

A.  BUOY  TENDING  OPERATIONS 

At  the  present  time,  buoy  tending  operations  are  at  best  an  inexact  science,  as  is  the 
case  in  many  shiphandling  routines.  The  manual  ship  positioning  techniques  of  experi¬ 
enced  deck  officers  have  been  raised  to  an  art  form.  Put  simply,  by  the  use  of  sextant 
angles  or  other  "local"  points  of  reference,  the  vessel  is  slowly  aimed  at  its  intended  po¬ 
sition.  As  the  distance  to  the  desired  location  decreases,  the  engines  are  slowed  or  * 
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Figure  1.  WLB  Comparison.  (From  Ref.  2) 


stopped  altogether.  As  the  vessel  crosses  this  desired  location,  the  buoy  is  released  in 
position  and  the  vessel  continues  on.  If,  however,  the  vessel  does  not  cross  the  desired 
buoy  position,  another  pass  must  be  made.  As  weather  conditions  deteriorate,  this 
procedure  can  rapidly  increase  on  station  time.  A  dynamically  positioned  vessel  can  be 
programmed  for  a  given  location  and  heading.  The  DP  vessel  adapts  to  environmental 
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changes  and  makes  optimal  use  of  thrusters  to  maintain  position  so  it  clearly  has  the 
advantage  ever  the  hit-or-miss  scenario  just  described. 
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111.  EQUATIONS  OF  MOTION 


The  mathematical  modeling  of  the  "plant",  the  waterborne  vessel  and  its  environ¬ 
ment,  is  divided  into  several  parts:  low  and  high  frequency  vessel  models,  wind  model, 
and  current  model. 

A.  COORDINATE  SYSTEM 

For  modelling  purposes,  two  different  coordinate  systems  are  used.  Earth-fixed  co¬ 
ordinate  axes  are  used  to  integrate  the  dynamic  equations  of  motion.  A  vessel  parallel 
set  of  axes,  with  the  same  origin  as  the  earth-fixed  axes,  is  used  for  computing  the  forces 
acting  on  the  vessel.  See  Figure  2  for  the  relationship  of  these  coordinate  axes  (Ref.  4: 
p.853).  The  x  and  y  axes  represent  the  vessel  in  surge  and  sway  and  the  angle  ^  is  the 
vessel  yaw  in  degrees  true.  Hence,  heading  angle  is  the  same  for  both  sets  of  axes. 


Figure  2.  Coordinate  Axes 
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B.  LOW  FREQUENCY  MODpL 

The  equations  of  motion  for  a  floating  vessel  in  a  seaway  are  highly  coupled  and 
non-linear.  Figure  3  is  a  block  diaqram  of  the  low  frequency  (LF)  motions  of  the  vessel 
in  surge,  sway,  and  yaw.  The  low  frequency  model  represents  the  dynamic  response  of 
the  ship  in  "calm"  water  with  wind,  current,  and  propulsive  force  disturbances.  The  ef¬ 
fect  of  wave  action  is  not  considered  here  as  that  is  modeled  by  the  high  frequency  model 
below.  The  coupled  non-linear  differential  equations  for  surge,  sway,  and  yaw  may  be 
expressed  as  follows: 


xL  =  1  lmx\my(S’L)iiL  +mxvc'l'L  -dx\xL-ue\ (xL -  uc)  +FWX  +Fp]  +  >iu  ( 1 ) 
h  =  1  lnty[myucilfL  -mx(xL  -uc)^L  -dy  |j;t  -  vc  |  {yL  -  vc)  +  Fwy  +Fr]  +  nLy  (2) 
4>l  =  l/”Vl  ”  H  ~  uc  Xfc.  —vc)  -d*\h\h  +W«  +MW  +  ^  (3) 


where 

•  *L  -  LF  surge  position  and  velocity 

•  3W’t  *  LF  sway  position  and  velocity 

•  'At.  i'L  •  LF  yaw  heading  angle  and  angle  rate 

•  we,  ve  -  Current  velocities  in  surge  and  sway 

•  Fw„  Fwy  -  Wind  force  in  surge  and  sway 

•  FP,  Ft  -  Prop  and  thruster  forces 

•  Mw,  Mt,  Me  -  Wind,  thruster,  and  current  moments  in  yaw 

•  i]Ly,  i}L*J/  -  Zero  mean  Gaussian  White  Noise  Processes 

•  d„  dy ,  d+  -  Coelficients  of  Drag 

•  m„  itty,  mt  -  Effective  Mass  Terms 

IRef.  4  :  p.  853] 
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Figure  3.  LF  Model 


C.  VESSEL  CHARACTERISTICS 

As  previously  stated,  the  modelled  vessel  is  a  USCG  WLB  "IRIS"  class.  The  vessel 
is  considered  to  have  the  following  characteristics: 

•  Length  (Overall)  -  180  feet 

•  Length  (between  perpendiculars)  -  170  feet 

•  Beam  -  37  feet 
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•  Draft  -  12  feet  \ 

•  Displacement  -  943.5  L.Tons 

•  lit,  =  71920  slugs 

•  nty  =  147615  slugs  '  ■ 

•  =  2.67x10'  slugfi1 

•  d,  =  169  slugjft 

•  dy  =  869  s/wg//f 

•  dt  -  1.17x10’  slugftjrad} 

[Ref.  3  :  p.  17] 

The  effective  mass  {m„  my,  m*)  terms  are  generated  by  a  computer  simulation  known 
as  the  Standard  Ship  Motion  Program  (SMP).  This  program  was  developed  at  the 
David  W.  Taylor  Naval  Ship  Research  and  Development  Center.  The  SMP  provides 
predictions  of  vessel  motion  in  translation  (surge  and  sway)  and  heading  (yaw)  in  irreg¬ 
ular  seas  based  upon  strip  theory.  The  effective  mass  is  the  sum  of  the  vessel  mass  (the 
mass  of  water  displaced  by  the  vessel)  and  an  "added  mass"  term.  The  added  mass  term 
is  predicted  by  the  SMP  versus  varying  frequencies  of  wave  encounter,  a.  Appendix  A 
is  the  printout  of  the  SMP  for  the  180  foot  WLB  buoy  tender.  The  added  mass  and  drag 
terms  shown  there  are  non-dimensional.  These  terms  are  a  function  of  the  frequency 
of  encounter  of  wave  motion,  a,  which  is  also  non-dimensional.  For  purposes  of  this 
thesis,  constant  dimensionalized  values  of  drag  and  mass  are  used,  based  upon  an  "av¬ 
erage"  wave  frequency  of  1.1  radians  per  second  which  corresponds  to  Sea  State  5  and 
an  average  wave  height  of  5.0  ft.  However,  the  drag  terms  of  the  SMP  cannot  be  di¬ 
rectly  related  to  the  LF  equations  of  motion  used  here.  The  drag  coefficients  are  derived 
directly  from  Bernoulli's  equations  for  fluid  flow  on  a  moving  object: 

F=  \l2pU1  (4) 

where  p  is  the  fluid  density  and  U  is  the  velocity  of  the  fluid  relative  to  the  object.  The 
drag  terms  are  the  product  of  the  Bernoulli  force  terms  with  the  effective  underwater 
area  of  the  vessel.  By  using  the  length  between  perpendiculars,  beam,  draft,  and  a  block 
coefficient  of  0.437,  the  drag  terms  are  as  listed  above.  Simulations  using  the  10DE 
software  package  confirm  these  values.  See  Appendix  B  for  the  simulation  results.  The 
vessel  is  simulated  at  zero  sea  state  with  maximum  thrust.  Using  the  effective  mass 
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terms  of  SMP,  the  drag  coelficients  are  then  varied  until  the  known  maximum  velocities 
arc  obtained. 


D.  HIGH  FREQUENCY  MODEL 

The  high  frequency  vessel  model  represents  the  high  frequency  wave  action  of  the 
vessel.  Three  harmonic  oscillators  are  modelled  in  surge,  sway,  and  yaw  with  variable 
frequency,  as  follows: 


XH-~  0>lxH  +  )JHx 

(5) 

ill  -  -  +  '///>■ 

(6) 

hi  -  ~ 

(7) 

d>x  =  '/WX 

(8) 

~  oy 

(9) 

~~  ^  w»#v 

(10) 

where 

•  x,„  x„  -III7  surge  position  and  velocity 

•  yW//  -  HE  sway  position  and  velocity 

•  ihh  'I'h  *  I  IF  yaw  angle  and  angle  rate 

•  «„  w ,,  wr  - 1  IF  angular  frequencies  in  surge,  sway,  and  yaw 

•  >L„  '/«,>»  '/«,.<■•  '/«,»  >lny>  *  Zero  mean  Gaussian  White  Noise  Processes 

|  Ref.  4  :  p.853] 

E.  WIND  MODEL 

The  wind  is  modelled  as  the  sum  of  slowly  varying,  low  mean  frequency  and  rapidly 
varying  high  frequency  gusts,  both  in  wind  speed  and  direction.  The  equations  for  the 
wind  model  arc: 


wvs  ~  '/  vm 


(ID 


**'i if  ~  a\  »V  "1”  '/ wy/ 


(12) 
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^ds  ~  *1wds 

(13) 

a2wdf+  Vwdf 

(14) 

|  where  the  model  terms  are  defined  as  follows: 

•  w„  -  Slowly  varying  wind  velocity 

•  wt/  -  Rapidly  varying  wind  velocity 

•  vv*  -  Slowly  varying  wind  direction 

•  vv u  -  Rapidly  varying  wind  direction 

•  >/*./.  >/»*>  -  Zero  mean  Gaussian  White  Noise  Processes 

The  wind  force  terms  of  equations  1 ,  2,  and  3  are  then  defined  by: 


!Twx  =  FMWvs  + 

(15) 

F\vy  ~  7*2 (/0(**Vf  wvf) 

(16) 

M*  =  W)Kf  +wv/)2 

(17) 

p  =  wds  +  Wdf-  (\pL  +  \jjn) 

(18) 

where 

•  rwy,  Mw  -  Wind  forces  and  moment 

•  p  -  Wind  angle  of  attack 

•  T7,,  Fj,  /j  -  Effective  wind  surface  areas 

[Ref.  4  :  p.  854] 

The  F,  terms  are  again  derived  from  Bernoulli's  equations  for  fluid  flow.  In  this  case 
the  fluid  is  air  and  the  area  is  the  effective  area  of  wind  resistance  of  the  vessel. 

F.  CURRENT  MODEL 

The  current  model  is  described  by  the  following  set  of  differential  equations: 


*C£  ~  *lxCE 

(19) 

S>CE  =  *1yCE 

(20) 
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(21) 


^C£  ”  *1*CE 


(22) 

(23) 

with  the  terms  defined  as  follows: 


uc  =  xCE  cos  \j/  +  yCE  sin  ^ 
ve  =  -*ce  sin  $  +ycE  cos  \J/ 


•  xcE>ycE  *  Current  velocities  in  earth  coordinates  (North  and  East  respectively) 

•  uc,  ve  -  Current  velocities  in  surge  and  sway 

[Ref.  4  :  p.  854] 

These  wind,  current,  high  frequency,  and  low  frequency  models  constitute  the  vessel 
"plant".  Each  model  represents  a  different  aspect  of  the  plant  dynamics.  It  is  clear, 
however,  that  the  only  model  which  may  be  controlled  by  the  DP  system  is  the  low  fre¬ 
quency  model. 


y 

i 
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IV.  CONTROL  SYSTEM  DESIGN 

A.  PROPULSORS 

The  design  of  the  buoy  tender  control  system  is  based  on  several  constraints.  First, 
there  are  only  three  propulsive  forces  available:  the  propeller  for  thrust  forward  and 
astern  (±or  direction,  surge)  and  the  bow  and  stern  thrusters  for  moment  control 
(±  ^direction,  yaw)  and  athwartships  motion  (±  y  direction,  sway). 

The  bow  and  stern  thrusters  are  200  HP  fixed  tunnel  units.  A  good  rule  of  thumb 
for  the  thrust  developed  by  fixed  azimuth  tunnel  thrusters  is  30  ft  IbJjSHP.  These  will 
yield  approximately  6000  pounds  force  in  the  ±  y  direction.  The  bow  thruster  is  located 
78.88  feet  forward  of  the  center  of  buoyancy  and  the  stern  thruster  is  81.12  feet  astern 
of  it.  See  Figure  4  for  the  thruster  layout. 

The  force  and  moment  equations  for  the  bow  and  stern  thrusters  are  thus: 


FT  =  DT+ST 

(24) 

MT  =  78.8857'-  81.12SF 

(25) 

where,, 

•  FT  -  Athwartships  force  due  to  thrusters 

•  MT  -  Yaw  moment  due  to  thrusters 

•  BT  -  Force  of  bow  thruster 

•  ST  -  Force  of  stern  thruster 

It  must  be  noted  that  athwartships  force,  FT,  is  assumed  to  be  applied  at  a  single  point, 
the  center  of  buoyancy.  The  moment  force,  MT,  is  applied  about  the  center  of 
buoyancy  by  the  thrusters  proportional  to  their  output  and  the  length  of  the  moment 
arm.  This  length  is  measured  from  the  center  of  buoyancy  to  the  thruster.  A  positive 
moment  is  considered  to  be  one  which  will  make  the  ship  turn  to  starboard,  a  negative 
moment  to  port.  See  Figure  4  for  the  thruster  moment  dimensions. 

The  fixed  pitch  propeller  is  driven  by  a  single  1170  SHP  diesel  electric  power  plant. 
It  can  normally  produce  upwards  of  20,000  pounds  force  of  thrust.  The  equation  for  the 
maximum  thrust  generated  by  the  propeller  is: 
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Figure  4.  Thruster  Forces  and  Moments 


FP  = 

1  1  max 


nmsnrmjiS(sr^f-\ 


max 


(26) 


Assume  a  propeller  efficiency,  tjf„p «  .75  and  a  maximum  velocity,  Kmix  =  14.3  knots ,  this 
yields  a  maximum  thrust  FPm„  =  20,250ft  Ibf.  However,  due  to  the  nature  of 
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positionkeeping,  the  propeller  controller  was  limited  to  12,000  pounds  force  in  the  DP 
mode. 

The  controller  design  is  based  only  on  the  low  frequency  vessel  model.  Responding 
to  the  high  frequency  wave  action  of  the  vessel  would  cause  the  controllers  to  modulate 
in  even  moderate  seas. 

The  low  frequency  system  was  placed  in  state  variable  format, 

x  =  Ax  +  Du  +  Gw  (27) 

y**Cx+  Du  +  v  (28) 

The  non-linear,  coupled  A  and  B  matrices  arc: 


0  —dx\x  —  uc\!mx  0  0 

0  0  0  1  0 

0  ~mx^lmy  0  -dy\y-  vc\lmy  0 

0  0  0  0  0 

0  -  {my  -  0  -(my-mx)xlm^  0 


0 

0 

0 

0 

0 

0 

1  l»ty 

1  ll»y 

0 

0 

78.88 

-81.12 

The  state  observations  matrix,  C  is: 

'l  0  0  0  0  0‘ 

C=  0  0  1  0  0  0 

0  0  0  0  1  0_ 

The  low  frequency  plant  is  linearized  about 

uc  —  vc  =  ijj  =  x  —  y  —  0 

and  the  vessel  parameters  substituted  to  yield  the  linear  A,  B,  and  G  matrices. 
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'0  1  0  0  0  0 

0  0  0  0  0  0 

0  0  0  1  0  0 

■  A  = 

0  0  0  0  0  0 

!  0  0  0  0  0  1 

i. 

0  0  0  0  0  0 


B 


0 

1.39£  — 05 
0 
0 
0 
0 


0 

0 

0 

6.77E-06 

0 

78.88 


0 

0 

0 

6.77£  — 06 
0 

-81.12 


The  plant  noise  matrix  is: 

0 

.0005 

0 

G  — 

.0005 

► 

0 

.0005 


The  state  variable  matrix  is: 

x 

x 

y 

X-  . 

y 

* 


(33) 


(34) 


(35) 


(36) 


« 


y 
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The  control  matrix  is:  , 

u  = 

The  Control  Toolbox  of  the  software  package  PC-MATLAB  (for  MS-DOS  Per¬ 
sonal  Computers,  Version  3.2-PC  June  8,  1987)  is  then  used  to  obtain  the  feedback 
gains.  Assuming  full  state  feedback  is  available,  the  lqr  (linear  quadratic  regulator) 
function  calculates  the  optimal  feedback  gain  matrix  K,p,  such  that  the  feedback  law: 


FP 

BT 

ST 


(37) 


minimizes  the  cost  function: 


J  = 


(xtQx  +  uT Ru)dt 


where  the  state  weighting  matrix,  Q  is: 


P 

Q  = 


3  £04  0 

0  3  £04 

0  0 

0  0 

0  0 

0  0 


0  0 

0  0 

6£04  0 

0  6  £05 

0  0 

0  0 


0  0 

0  0 

0  0 

0  0 

1  £09  0 

0  1  £09 


and  the  control  weighting  matrix,  R  is: 


R  = 


1  0  0 
0  1  0 
0  0  1 


(39) 


(«) 


(41) 


» 
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The  lqr  function  returns  the  follqwing  feedback  gains: 


173  4977  0  0  0  0 

.0  0  150  4434  24962  249962 

'0  0  193  5701  -19415  -19415 


This  places  the  eigenvalues  of  the  closed  loop  linearized  system  at: 

" -0.035  ± 7'0.035 ' 

—0.034  ±70.034 

-1.0  (43) 

-3.6£06 

Since  all  the  states  are  not  available  for  feedback,  an  observer  is  utilized  to  provide 
an  estimate  of  the  states.  A  Kalman  filter  is  designed  with  the  aid  of  the  PC-MATLAB 
function  lqe  (linear  quadratic  estimator).  Recall  the  system: 

x  =  Ax  +  Bu  +  Gw  (44) 


y  =  Cx  +  Du  +  v  (45) 

* 

where  the  expected  values  of  the  plant  noise,  w,  and  the  measurement  noise,  v,  are 
£[u»]  =  £[v]  =  0  and  the  covariance  matrices  for  plant  and  measurement  noise  are,  re¬ 
spectively,  £[wu>']  and  ££w'].  The  covariance  matrices  are: 

£[vw']  =  0.0005  (46) 


'100  0  0 
E[_wy  —  0  100  0 

0  0  .0003 


The  lqe  function  returns  the  Kalman  filter  gains,  K,  for  the  equation: 


x  =  Ax  +  Bu  +  K(y  —  Cx) 


where  x  is  an  estimate  of  x.  The  Kalman  gain  matrix,  K,  is: 
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0.0669  0  0 

0.0022  0  0 

0  0.0669  0 

K  = 

0  0.0022  0 

0  0  1.6069 

0  0  1.2910 

The  output  of  the  Kalman  filter  is  the  sum  of  the  nonlinear  low  frequency  system,  in¬ 
cluding  sampled  values  of  position,  heading,  and  wind  velocity  and  direction,  with  the 
Kalman  updates  of  position.  The  sampled  wind  velocity  and  direction  data  is  trans¬ 
formed  into  ship's  coordinates  and  forces  the  vessel  Kalman  filter  exactly,  within  sam¬ 
pling  error,  as  the  true  wind  forces  the  low  frequency  model.  The  optimal  feedback  gains 
are  then  used  with  the  Kalman  filter  output  to  create  the  optimal  Kalman  filtered 
control: 


(49) 


»/  A 

u  — 

The  reference  error  matrix,  e,  is  defined  by: 

A 

rx~x 

A 

ri~X 

A 

a  ry-y 

<?  =  A 

ry-y 

A 

The  r  values  here  represent  the  reference  positions,  heading,  and  associated  velocities. 
The  reference  positions  and  heading  are  the  desired  location  of  the  vessel.  Since  the 
desired  location  is  to  be  maintained,  clearly  the  reference  velocities  must  be  zero. 

These  optimal  controls  are  then  passed  to  a  non-linear  filter.  This  non-linearity  is 
a  limiter  which  prevents  the  thrusters  and  main  propulsion  from  exceeding  the  maximum 
DP  mode  limits  (i.e.,  12,000  Ibft or  main  propulsion  and  6,000  /6/for  thrusters ). 

It  is  readily  apparent  that  the  control  system  as  designed  does  not  correct  for  the 
steady  state  error  caused  by  the  introduction  of  current.  A  modification  to  the  control 
design  to  compensate  for  current  offset  is  imperative.  A  Kalman  filter  predictor  of 


(50) 


(51) 
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current  is  designed  based  upon  the  residual  error  in  the  vessel  position  estimator.  That 
is,  the  difference  in  estimated  and  measured  position  and  heading  is  assumed  to  be  due 
to  current.  This  is  the  "measured"  value  of  the  Kalman  current  filter.  Recall  the  Kalman 
filter  equation: 


x  —  Ax  +  Du  +  tf(y  CJr) 


For  the  current  predictor  filter: 


A 

X  ?= 


A 


A 


A 

Mc 


y  = 


0  1  0  0  0  0 
0  0  0  0  0  0 
0  0  0  1  0  0 
0  0  0  0  0  0 
0  0  0  0  0  1 
0  0  0  0  0  0 


A 

xT  —  X 

A 

yr-y 

A 

\j/r- i> 


c= 


1 

0 

0 


0  0 
0  1 
0  0 


0  0  0 
0  0  0 
0  1  0 


(52) 


(53) 


(54) 


(55) 


(56) 
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The  current  estimator  covariance  piatrices  arc: 


E\_ww,m]c  —  0.00002 


t 

£[vv']c 


"5  0  0 

0  5  0 

0  0  .005 


The  current  estimation  filter  gain  matrix,  Ktt  is: 


Kc  = 


0.0140 

0.0044 

0 

0 

0 

0 


0 

0 

0.0140 

0.0044 

0 

0 


0 

0 

0 

0 

6.0000 

18.200 


(57) 

.(58) 


(59) 


The  terms  xT,yT,  and  are  the  observations  of  surge,  sway,  and  yaw  including  high 
frequency  wave  action  and  measurement  noise.  The  estimates  of  current  are  converted 
to  force  terms  by  Bernoulli's  equation  and  subtracted  from  the  appropriate  thruster 
control  commands.  Figure  5  shows  the  final  system  control  design  in  block  diagram 
form.  ? 
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V.  COMPUTER  SIMULATION 


The  equations  of  motion  and  control  system  discussed  Hn  previous  chapters  tvere 
implemented  using  the  Dynamic  Simulation  Language/ VS  (DSL).  It  is  a  high  level 
continuous  simulation  language  developed  by  IBM  which  incorporates  VS  FORTRAN 
as  a  subset.  The  program  is  included  in  Appendix  C  here  as  an  immediate  reference. 
Comments  are  included  as  additional  documentation.  DSL  was  used  here  since  it  is  a 
high  level  language  which  models  the  system  in  continuous  time.  Although  the  system 
is  modeled  as  an  analog  system,  DSL  allows  for  the  sampling  of  data  as  well.  This 
sampling  was  used  with  the  position,  heading,  and  wind  measurements,  which  were  up- 
dated  at  1  second  intervals.  Nearly  all  of  the  various  models  require  Gaussian  noise  in¬ 
puts  which  arc  represented  by  the  NORMAL  function.  This  function  has  inputs  of  seed 
number,  mean,  and  standard  deviation.  One  final  word  on  DSL  is  that  it  has  its  own 
plotting  routine.  System  outputs  presented  here  were  printed  using  this  plotting 
package. 
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VI.  RESULTS 

The  simulation  results  of  the  dynamically  positioned  buoy  tender  are  shown  here  in 
a  number  of  different  environments.  As  discussed  previously,  the  specifications  for  the 
desired  dynamic  positioning  system  arc  ambitious.  To  review,  the  vessel  is  to  attain  and 
maintain  position  to  within  a  10  meter  radius  in  30  knots  of  wind  in  any  aspect  and  5 
knots  of  current  within  10  degrees  of  the  bow.  The  initial  test  of  the  simulation  is  with 
no  current,  no  wind,  and  at  sea  state  0.  These  conditions  introduce  the  least  amount 
of  disturbance  into  the  system.  As  can  be  seen  from  Figure  6,  the  vessel  holds  position 
in  surge  to  within  5  feet.  Similarly,  Figure  7  and  Figure  8  show  the  vessel  holding  sway 
position  o  within  5  feet  and  a  decaying  exponential  response  in  yaw  angle,  with  its 
maximum  overshoot  of  0.05  radians.  In  each  graphical  presentation  of  vessel  position 
and  heading,  the  true-measured  position  and  heading  are  plotted  along  with  the  respec¬ 
tive  vessel  Kalman  filter  estimates.  A  position  plot  is  shown  in  Figure  9.  The  position 
plot  is  presented  as  an  "overhead”  look  at  the  vessel  track  in  vessel  coordinates,  in  order 
to  obtain  this  vessel  response,  a  conservative  use  of  propulsive  force  is  used,  as  witnessed 
in  Figure  10  and  Figure  11. 

The  system  is  again  simulated  at  sea  state  0  with  no  wind  and  a  current  of  3  feet  per 
second/  In  Figure  12  the  system  performance  is  observed  in  the  surge  direction.  Since 
the  vessel  estimator  is  not  forced  by  the  current  estimator,  there  is  a  bias  introduced  by 
the  vessel  estimator.  This  bias  could  be  corrected  by  modifying  the  state  variable  matrix 
with  the  current  velocities  in  surge  and  sway  such  that: 


This  bias  is  exactly  balanced  in  the  propulsion  scheme  by  the  current  estimated  force. 
Similarly,  the  system  performance  in  sway,  Figure  13,  and  yaw,  Figure  14,  is  within  the 
system  constraints  in  the  presence  of  this  current.  The  position  plot,  Figure  15,  closely 
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resembles  its  counterpart  for  the  zero  current  ease.  However,  the  overshoot  in  both 
surge  and  sway  is  even  more  apparent  in  this  plot.  The  transient  performance  is  not  of 
particular  importance  us  the  steady  state  performance  is  acceptable.  Due  to  the  con¬ 
tinuing  nature  of  the  current  disturbance,  non-zero  propulsive  force.,  result  in  the  steady 
state  as  seen  in  figure  16  for  main  propulsion  force  and  in  Figure  i7  for  bow  and  stern 
thruster  forces. 

The  vessel  when  simulated  at  a  current  velocity  of  4  feet  per  second  docs  not  main¬ 
tain  its  desired  position  due  to  the  steady  state  saturation  of  the  bow  and  stern  thrusters. 
The  maximum  current  for  which  this  system  is  effective  is  considered  to  be  3  feet  per 
second. 

Since  fully  developed  seas  arc  correlated  with  wind  velocity,  the  final  simulation  is 
at  sea  state  5  with  current  of  3  feet  per  second  and  30  knot  winds.  With  the  addition 
of  the  Kalman  predictor  for  current  and  the  fact  that  wind  is  sampled  each  second  and 
is  incorporated  into  the  vessel  position  Kalman  filter,  no  great  degradation  in  system 
performance  is  anticipated.  The  surge,  sway,  and  yaw  performance  of  the  vessel  can  be 
seen  in  Figure  18,  Figure  19,  and  Figure  20,  respectively.  The  vessel  response  in  surge 
and  sway  are  well  within  the  specifications.  Although  thers  is  no  requirement  to  main¬ 
tain  vessel  heading,  the  steady  state  yaw  response  of  the  vessel  is  held  to  about  0.2  radian 
of  the  reference.  The  load  on  the  propulsors  under  these  conditions  is  heavy  as  can  be 
seen  in  Figure  21  and  Figure  22.  The  position  plot  of  the  vessel  under  these  conditions 
is  seen  in  Figure  23.  The  steady  state  position  is  still  within  the  vessel  specifications. 

Although  there  is  a  degradation  in  system  performance  with  increased  environ¬ 
mental  disturbances,  the  vessel  is  able  to  maintain  a  desired  position  to  within  the  con¬ 
straints  imposed  at  the  outset. 
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Figure  6.  Response  in  Surge:  Sea  state  0,  No  current,  No  wind.  XT  is  measured 
surge  position  and  XE  is  estimated  surge  position. 
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Figure  7.  Response  in  Sway:  Sea  state  0,  No  current,  No  wind.  YT  is  measured 
sway  position  and  Y’E  is  estimated  sway  position. 
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Figure  8.  Resnoiv  ?  in  Yaw;  Sea  state  0,  No  current,  No  wind.  PT  is  measured 
yaw  angle  and  PE  is  estimated  yaw  angle. 
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Figure  12.  Surge  Response:  Sea  state  0,  No  wind,  Current  =  3  fps.  XT  is 
measured  surge  position  and  XE  is  estimated  surge  position. 
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Figure  13.  Sway  Response:  Sea  state  0,  No  wind,  Current  =  3  fps.  YT  is 
measured  sway  position  and  YE  is  estimated  sway  position 

y 
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Figure  14.  Yaw  Response:  Sea  state  0,  No  wind,  Current  =  3  fps.  PI  is  meas¬ 
ured  yaw  angle  and  PE  is  estimated  yaw  angle. 
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Figure  15.  Position  Plot:  Sea  state  0,  No  wind,  Current  =  3  fps 
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Figure  16.  Main  Propulsion  Force:  Sea  state  0,  No  wind,  Current  =  3  fps 
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Figure  17.  Bow  and  Stern  Thruster  Forces :  Sea  state  O,  No  wind,  Current  =  3 
fps 
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Figure  19.  Sway  Response:  Sea  state  5,  Wind  =  30  knots,  Current  =  3  fps.  YT 
is  measured  sway  position  and  Y'E  is  estimated  sway  position. 
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Main  Propulsion  Force:  Sea  state  5,  Wind  =  30  knots,  Current  =  3 
ips 
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Figure  22.  Bow  and  Stern  Thruster  Forces:  Sea  state  5,  Wind  =  30  knots,  Cur¬ 
rent  =  3  fps 
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VII.  CONCLUSIONS 


This  computer-based  dynamic  positioning  system  is  a  viable  alternative  to  manual 
position  keeping  methods  in  environmental  conditions  of  up  to  sea  state  5.  current  up 
to  1.5  knots,  and  wind  up  to  30  knots,  particularly  when  precision  requirements  arc 
deemed  to  be  of  greater  importance  than  fuel  cost.  In  order  to  meet  the  requirments  for 
the  nev  buoy  tender  licet,  bow  and  stern  thrusters  of  400  SUP  should  be  installed. 
These  will  allow  the  DP  system  to  maintain  vessel  position  and  heading  in  the  environ¬ 
ments  specified.  In  addition,  variable  azimuth  thrusters  will  provide  enhanced  control 
system  performance  over  fixed  tunnel  thrusters. 
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APPENDIX  A.  STANDARD  SHIP  MOTION  PROGRAM 


The  added  mass  and  drag  coefficients  for  the  ISO  foot  WLB  "IRIS"  class  buoy 
tender  arc  dimensionalized  forms  ol  the  values  presented  here. 


II*  IhI  Nt*  I’WICI  l*l»  ClAt**l  I*  J/A*  |*lt 
HKt  or  INI*  »A*IKVU»I 


INI*  CHiUCttlllMCI  * 

|N|*  UNIIN  (Iff) 
Ml*  *1  NtMNl** 
HIM  At  N|MH|*| 

MtrtACtMil  ll.tr.l 
it  I ion  »H|*  «*((* 


VlltlCH  l*t*f I0»l  * 

e.  or  imvim  tvcoi*  t.ir  'in 

c.  it  omviti  ircim  m. 77  mi 

NCUttMl*IC  Hi,  l*NI  i.tr  i  it  i 

NttACMHt*  11*111  17.14  lltl 

c.  or  ivoiahct  inti*  7.i*  rttt 


LO«lttvMO*t  lOCAIIOH***  * 

e.  or  tiivnr  <ttoi  *>.*»  rttt 
c.  or  ivHiKi  utii  inn  rttt 
c.  or  neiAttOM  tun  •*.**  rttt 


NO  1 1  IN  CH***C1tMITlt»  • 

•oil  otuotut  t*.»*  rttt 

r  I icm  oiiamui  <7.5*  rtri 

tam  *i**»iv«  ii. s*  rttt 

tsiiNiit*  »ou  rtnor  io.li  ttcom* 


to*ruiti  **t  a*  - 

VAIttriAMt  4J4I.7  10.  lltl  Mf/(l****C*M»  ♦.*»* 

btmt  tuiriett  nuu  iuj.i  to.  rttt  vi/oimzimih  ».W 


nuu  cotmcttHii  • 

not*  KM  M)> 

ItCMO*  ttll  *.775 

mtNAtlC  K*l  *.5*4 


«  yiimiM  AtrtitMtt 

•t  itti  itrtrtNtt 
•  Mf,7.  ItrtttNCt 


ICO/ltNUN  *.  4t1 
IC>/ltN6lH  *.*») 
Icr/ifNOiH  *.56* 


RI/ltAN  *.46* 

ri/u*  • .  7!* 

t*/i*  *  *.::» 

(Oil  HtO  t|A*|AN|t  «.S*I 


VCC/M»n  *.**! 

M/*tAN  *  • • J*4 

ON/HAN  *.071 

*  A/M  AN  *.*4) 

•  l/lf  AN  *.l*4 


17*.**  Ittl  ItNIlH/ll AN  4.5*5 

)7.l*  rttt  It AN/ |AAft  5. 1*1 

ti.**  rttt  ham/man  *.s:t 

HJ.5l.tONt  *11*1/1.  *11**14*5  1*7. *4* 

ii.  •*  *  no  i*  reottre  NVNit*  *.:7i 


it*  r«*t  wn  cuica  i*ii  ctAitn  i*  j**. •  i*i* 

t All t  0*  $N|*  AirtNlAOt  rAAtlCMAtt 


i*to  tNA*Atiminti  - 

IMO  ItNOIH  AlONI  lltl  Ittl  NO.  ||  11.50  rttl 

lltl  M| | QHl  (to  *0.  II  16.61  rttl 

ioiai  vim*  te*r Act  A*t*  i*«to  itt  no.  it  jia  in.  rttt 

A*|A  /  Itn  I  MAfll  0.077 


Aunt*  CNAiActmitic*  • 

•Ullt *  «0«f  t HO* I  UNOtH  mi  no.  \t  1.75  rttt 

lUflt*  t|l  CNCtl  Ltnoib*  lltl  NO,  I)  1.75  Met 

•gilt*  *t an  l* an  (Ml  NO.  II  1,25  rttt 

i*iai  **t  1 1 1  *  igifAtt  Ait  a  Hunt*  itt  *o,  ii  J7  io.  rttt 

.  »  Alt  A  /  U**  I  MAM)  *.**< 


nok i  ir  a  'itt*  ttntitNti  a  r a i *  or  Ammoti  *i.i.»  iu*r  ittin. 
twtN  IMt  Itllll  *l7NrAtt  II  Confute*  r*»  fwt  ioiai  ama  or  ioin 


Figure  24.  SSMP  Vessel  Particulars 
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Figure  25.  SSMP  Added  Mass  and  Damping  Coefficients 
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APPENDIX  B.  IODE  SIMULATIONS 


The  following  are  printouts  from  IODE  programs.  IODE  (Interactive  Ordinary 
Differential  Equations)  is  an  interactive  software  package  that  runs  on  the  VM/CMS 
mainframe.  IODE  interrogates  the  user  about  the  system  of  dillcrcntial  equations  and 
provides  the  solution  in  tabular  and/or  graphical  form. 

The  first  program  of  this  appendix  is  a  model  of  the  low  frequency  vessel  in  surge. 
The  vessel  is  modeled  at  full  propulsive  force,  with  no  environmental  disturbances.  Us¬ 
ing  the  SMI’  value  for  vessel  effective  mass  and  the  derived  value  for  the  drag  coefficient, 
the  vessel  is  simulated  in  the  open  loop  configuration  and  velocity  plotted  against  time 
in  figure  26  on  page  50  The  IODE  program  listing  is: 

VARIABLES  &  INITIAL  CONDITIONS: 

X  =  .  0 
XD  =  .  0 
T  =  .  0 

CONSTANTS: 

MX  =  71920. 00000 
NX  =  .0 

DX  =  169. 0000000 
FWX  =  . 0 
UC  =  .  0 

U  =  20250. 00000 

SPECIAL  FUNCTIONS: 

CC  =  XD-UC 

DERIVATIVES: 

D(X  /D(T  )  =  = 

XD 

D(XD  /D(T  )  =  = 

NX+( U+FWX-DX*ABS( CC)*CC ) /MX 

OUTPUTS: 

TITLE:  SURGE  NO  FEEDBACK 
NO  TABULATION 
PLOT:  X  XD 

AGAINST:  T 
PLOT:  XD 

AGAINST:  T 


AT  INTERVAL  .2000000000 


*  XD 


10.9  4 


7.29 


5.47 


0.  OOOE+OO 
-4 - 4- 


57.5 

-4 - 4- 


75.0 

-4- - 4- 


115. 

-4  —  - 


9.11  4 


5.65 


1.82 


0 . 000E+004  . 

--4------- 

T  0 . OOOE+OO 


37.5 

SURGE  NO  FEEDBACK 
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75.0 


-4 - 

115. 


150. 


4  10.9 
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4  0. OOOE+OO 
-4— 

150. 


Figure  26.  Surge  Motion  with  No  Feedback 
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As  in  the  first  simulation,  thp  low  frequency  sway  of  the  vessel  is  modeled  in  the 
open  loop  and  velocity  plotted  against  time  in  Figure  27  on  page  51.  The  IODE  pro¬ 
gram  for  sway  is: 


VARIABLES  &  INITIAL  CONDITIONS: 
Y  =  .0 
YD  =  .0 
T  =  .  0 

CONSTANTS: 

MY  =  147615.  0000 
DY  =  869. 0000000 
FWY  =  .  0 
VC  =  .0 

FT  =  12000. 00000 
NY  =  .  0 


SPECIAL  FUNCTIONS: 

CC  =  YD-VC 

DERIVATIVES: 

D(Y  /D(T  )  =  = 

YD 

D(YD  /D(T  )  =  = 

NY+((FT+FWY-DY*ABS(CC)*CC)/MY) 

OUTPUTS: 

TITLE:  SWAY  NO  FEEDBACK 
NO  TABULATION 
PLOT:  YD 

^AGAINST:  T  AT  INTERVAL  .2000000000 
END  CALCULATION  WHEN  T  .  GE.  150. 000 
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150. 
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Figure  27.  Sway  Motion  with  No  Feedback 


Finally,  the  yaw  motion  of  the  low  frequency  model  is  simulated  in  the  open  loop 
mode  and  heading  angle  rate  plotted  against  time  in  Figure  28. 
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VARIABLES  &  INITIAL  CONDITIONS: 

P  =  .  0 
PD  =  .  0 
T  =  .  0 

CONSTANTS: 

MP  *  267000000.0 
DP  =  1170000000. 

NP  =  .0 
MW  =  .0 
MC  =  .  0 

MT  =  954000.  0000 
DUM  =  .  0 

DERIVATIVES: 

D(P  /D(T  )  =  = 

PD 

D(PD  /D(T  )  =  = 

NP+( MT+MW+MC-DP*ABS( PD)*PD) /MP 

OUTPUTS: 

TITLE:  PS1  DOT 
NO  TABULATION 
PLOT:  PD 

AGAINST:  T  AT  INTERVAL  .2000000000 
END  CALCULATION  WHEN  T  .  GE.  150. 000 


C 
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Figure  28.  Yaw  Motion  with  No  Feedback 
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APPENDIX  C.  PROGRAM  LISTING 

The  following  is  a  listing  of  the  DSL  program  which  simulates  the  vessel,  its  envi¬ 
ronment,  and  the  dynamic  positioning  control  system.  Comments  arc  included 
throughout  to  augment  the  explanations  provided  in  the  body  of  the  thesis. 


’•'DEVELOPED  BY  W.  R.  CAIRNS,  LCDR,  USCG  FOR  THESIS  RESEARCH: 

*"DESIGN  AND  SIMULATION  OF  A  DYNAMIC  POSITIONING  SYSTEM  FOR  A  U.  S. 
’•'COAST  GUARD  BUOY  TENDER" 

TITLE:  DPSIM5 

’’'THIS  PROGRAM  SIMULATES  A  WLB  "IRIS"  CLASS  USCG  BUOY  TENDER  SUBJECT  TO 
’•'WIND,  CURRENT,  AND  WAVE  FORCES  AND  CONTROLLED  BY  A  DYNAMIC  POSITIONING 
’•'SYSTEM 

*THE  PROGRAM  PARAMETERS  (PARAM)  ARE  DEFINED  AS: 

*  FI ,F2 ,F3  ARE  THE  VESSEL  EFFECTIVE  WIND  AREAS  IN  SURGE, SWAY,  AND  YAW 

*  A1,A2  ARE  THE  VESSEL  WIND  VELOCITY  AND  DIRECTION  FORCE  COEFFICIENTS 

*  RX,RY,RP  ARE  THE  REFERENCE  POSITIONS  IN  SURGE, SWAY  AND  YAW 
PARAM  Fl=2. 294,F2=11. 16,F3=474.  3,. .  . 

A1=0. 001 ,A2=0.  001,. . . 

W=00.  8,.  . . 

IWVS=30.  0,.  .  . 

IWDS=0.  0,.  .  . 

XCI=0.  0,.  . . 

YCI=0. 0,. . . 

PCI=0. 0,. . . 

RX=100.  0,.  .  . 

RY=100. 0, . . . 

RP=1. 00 

*THE  CONSTANTS  (CONST)  ARE  DEFINED: 

*  X0,XD0  INITIAL  CONDITIONS  FOR  SURGE  VELOCITY  AND  ACCELERATION 

*  Y0,YD0  INITIAL  CONDITIONS  FOR  SWAY  VELOCITY  AND  ACCELERATION 

*  PO.PDO  INITIAL  CONDITIONS  FOR  YAW  VELOCITY  AND  ACCELERATION 

*  MX , DX  MASS  AND  DRAG  COEFFICIENTS  IN  SURGE 

*  MY , DY  MASS  AND  DRAG  COEFFICIENTS  IN  SWAY 

*  MP,DP  MASS  AND  DRAG  COEFFICIENTS  IN  YAW 

*  NX , UX , SX  SEED, MEAN, STD  DEVIATION  FOR  LF  PLANT  NOISE  IN  SURGE 

*  NY,UY,SY  SEED, MEAN, STD  DEVIATION  FOR  LF  PLANT  NOISE  IN  SWAY 

*  NP,UP,SP  SEED, MEAN, STD  DEVIATION  FOR  LF  PLANT  NOISE  IN  YAW 

*  HNX,HUX,HSX  SEED, MEAN, STD  DEVIATION 

*  HNY,HUY,HSY  SEED , MEAN , STD  DEVIATION 

*  HNP,HUP,HSP  SEED, MEAN, STD  DEVIATION 

*  NWX , UWX , SWX  SEED, MEAN, STD  DEVIATION 

*  NWY , UWY , SWY  SEED, MEAN, STD  DEVIATION 

*  NWP,UWP,SWP  SEED, MEAN, STD  DEVIATION 

*  MNX , UMX , SMX  SEED, MEAN, STD  DEVIATION 

*  MNY , UMY , SMY  SEED, ME AN, STD  DEVIATION 

*  MNP,UMP,SMP  SEED, MEAN, STD  DEVIATION 

*  CNX,CUX,CSX  SEED, MEAN, STD  DEVIATION 

*  CNY,CUY,CSY  SEED, MEAN, STD  DEVIATION 

*  CNP,CUP,CSP  SEED, MEAN, STD  DEVIATION 

*  NCX,MCX, SCX  SEED, MEAN, STD  DEVIATION 


FOR  HF  PLANT  NOISE  IN  SURGE 

FOR  HF  PLANT  NOISE  IN  SWAY 

FOR  HF  PLANT  NOISE  IN  YAW 

FOR  HF  FREQ  NOISE  IN  SURGE 

FOR  HF  FREQ  NOISE  IN  SWAY 

FOR  HF  FREQ  NOISE  IN  YAW 

FOR  LF  MEAS  NOISE  IN  SURGE 

FOR  LF  MEAS  NOISE  IN  SWAY 

FOR  LF  MEAS  NOISE  IN  YAW 

FOR  CURRENT  NOISE  IN  SURGE 

FOR  CURRENT  NOISE  IN  SWAY 

FOR  CURRENT  NOISE  IN  YAW 

FOR  CURRENT  NOISE  IN  SURGE(KALMAN) 
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*  NCY,MCY,SCY  SEED, MEAN, STD  DEVIATION  FOR  CURRENT  NOISE  IN  SWAY  (KALMAN) 

*  NCP,MCP,SCP  SEED, MEAN, STD  DEVIATION  FOR  CURRENT  NOISE  IN  YAW  (KALMAN) 

*  NWVS,UWVS,SWVS  SEED, MEAN, SD  FOR  SLOW  WIND  VELOCITY  NOISE 

*  NWVF,UWVF,SWVF  SEED, MEAN, SD  FOR  FAST  WIND  VELOCITY  NOISE 

*  NWDS , UWDS , SWDS  SEED, MEAN, SD  FOR  SLOW  WIND  DIRECTION  NOISE 

*  NWDF,UVDF,SWDF  SEED, MEAN, SD  FOR  FAST  WIND  DIRECTION  NOISE 

*  FPR,FPL  PROPULSION  LIMITER  (+/-) 

*  BTR , BTL  BOW  THRUSTER  LIMITER  (+/-) 

*  STR , STL  STERN  THRUSTER  LIMITER  (+/") 

*  EXD,EXU  DEAD  SPACE  IN  SURGE 

*  EYD,EYU  DEAD  SPACE  IN  SWAY 

*  EPD,EPU  DEAD  SPACE  IN  YAW 

*  RXD,RYD,RPD  REFERENCE  VELOCITIES  IN  SURGE,  SWAY,  YAW 

*  K  VESSEL  KALMAN  FILTER  GAINS  SUBSCRIPTED  ACCORDING  TO  MATRIX 

*  C  CURRENT  ESTIMATION  FILTER  GAINS  SUBSCRIPTED  ACCORDING  TO  MATRIX 

*  G  OPTIMAL  FEEDBACK  GAINS  SUBSCRIPTED  ACCORDING  TO  MATRIX 
CONST  X0=0. 0,XD0=0.  0,. .  . 

Y0=0. 0, YD0=0.  0,.  .  . 

P0=0.  0 , PD0=0. 0,. . . 

MX=71920.  0,DX=169. 00000,. . . 

MY=147615. 0,DY=869. 0000,. . . 

MP=2.  67D+08 ,DP=1. 17D+09,. . . 

NX=19.  0 ,UX=0.  0,SX=0. 0005,. . . 

NY=2 1.0,  UY=0.  0 ,  SY=0.  0005 ,..  . 

NP=23.  0,UP=0.  0,SP=0. 00001,. . . 

HNX=19.  0,HUX=0. 0,HSX=3. 0,. . . 

HNY=21.  0,HUY=0.  0,HSY=3. 0,. . . 

HNP=23.  0 ,11UP=0.  0,HSP=0.  175,.  .  . 

NWX=19.  0,UWX=0.  0,SWX=0.  005,.  .  . 

NWY=21. 0,UWY=0. 0,SWY=0. 005,. . . 

NWP=23. 0,UWP=0. 0,SWP=0. 000175,. . . 

MNX=19. 0,UMX=0. 0,SMX=3. 0000,. . . 

MNY=21.  0,UMY=0.  0,SMY=3. 0000,. . . 

MNP=23.  0,U!1P=0.  0,SMP=0.  0175,.  .  . 

CNX=19. 0,CUX=0. 0,CSX=0. 00050,. . . 

CNY=21. 0,CUY=0. 0 ,CSY=0. 00050,. . . 

CNP=23.  0,CUP=0. 0,CSP=0. 00001,. . . 

NCX=19. 0,MCX=0. 0,SCX=0. 0005,. . . 

NCY=21. 0,MCY=0. 0,SCY=0. 0005,. . . 

NCP=23. 0 ,MCP=0. 0,SCP=0. 00001,. . . 

NWVS=19.  0,UWVS=0. 0,SWVS=0.  0005,.  .  . 

NWVF=2 1 . 0,UWVF=0. 0,SWVF=0. 0005,. . . 

N\v’DS=23.  0,UWDS=0.  0,SWDS=0.  00001,.  .  . 

NWDF=23.  0,UWDF=0.  0,SWDF=0.  00001,.  .  . 

VNN=15. 0, VNU=0. 0, VNS=0. 0005,. . . 

FPR=12000. 0,FPL=-12000. 0,. . . 

BTR=6000.  0,BTL=-6000.  0,. .  . 

STR=6000. 0,STL=-6000. 0,. . . 

RXD=0. 0,RYD=0. 0,RPD=0. 0,.  .  . 

K11=0. 0669 ,K12=0. 0,K13=0. 0,. . . 

K21=0. 0022,K22=0. 0,K23=0. 0,. . . 

K31=0. 0,K32=0. 0669,K33=0. 0,. .  . 

K41=0. 0 ,K42=0. 0022,K43=0. 0,. .  . 

K51=0. 0,K52=0. 0,K53=1. 6069,. . . 

K61=0. 0,K62=0. 0,K63=1.  2910,.  .  . 

C11=0. 014,C12=0.  0,C13=0.  0,.  .  . 
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C21=0. 0044,C22=0. 0,C23=0.  0,.  .  . 

C31=0. 0,C32=0.  014,C33=0.  0,.  .  . 

C41=0. 0,C42=0. 0044,C43=0.  0,.  .  . 

C51=0. 0,C52=0. 0,C53=6.  0000,.  .  . 

C61=0. 0 ,C62=0.  0,C63=18.  200,.  .  . 

Gll=193. 0,G12=5701. 00,G13=0. 0,G14=0. 0,G15=0. 0,G16=0. 0,. . . 

G21=0. 0,G22=0. 0,G23=150.  0,G24=4434.  0,G25=24962.  0,G26=24962.  0,.  . . 
G31=0. 0,G32=0. 0,G33=193.  0,G34=5701.  0,G35=-19415.  000 ,G36=-19415.  00 

*  Gll=78. 50,G12=9200. 00,G13=0. 0,G14=0. 0,G15=0. 0,G16=0. 0,. . . 

*  G21=0. 0,G22=0. 0,G23=42. 00,624=7200.  0,G25=5854.  0 ,G26=8200.  0 , .  .  . 

*  G31=0. 0 ,G32=0. 0,G33=38. 0,G34=6900.  0,G35=-4553.  000,G36=-7S00.  00 

*  Gll=77. 46,G12=3335. 50,G13=0.  0,G14=0.  0,G15=0.  0,G16=0.  0,.  . . 

*  G21=0.  0 ,G22=0. 0,G23=47. 556,G24=2483. 80,G25=249. 62,G26=136. 7,. . . 

*  G31=0. 0,G32=0. 0,G33=61.  14,G34=3193.  40,G35=-194.  150,G36=-106.  34 
DYNAMIC 

,vL0W  FREQUENCY  MODEL  PLANT  NOISE  IN  SURGE,  SWAY,  AND  YAW 
NOIX=NORMAL( NX , UX , SX) 

NOI  Y=NORMAI,(  NY ,  UY ,  SY) 

NOIP=NORMAL( NP , UP , SP) 

*HIGH  FREQUENCY  MODEL  PLANT  NOISE  IN  SURGE,  SWAY,  AND  YAW 
HNOIX=NORMAL( HNX , HUX , USX) 

HNOI Y=NORMAL( HN Y , HUY , HS Y ) 

UNO I P=NORMAL( HNP,HUP,HSP) 

’•'MEASUREMENT  NOISE  IN  SURGE,  SWAY,  AND  YAW 
MNOIX=NORMAL( MNX , UMX , SMX ) 

MNOI Y=NORMAL( MNY , UMY , SMY) 

MNOI P=NORMAL( MNP , UMP , SMP) 

*CURRENT  NOISE  IN  SURGE,  SWAY,  AND  YAW 
CXNOI=NORMAL( CNX , CUX , CSX ) 

CYNOI=NORMAL( CNY , CUY , CSY) 

CPNOI=NORMAL( CNP , CUP , CSP) 

*KALMAN  FILTER  CURRENT  NOISE  IN  SURGE,  SWAY,  AND  YAW 
KCX=NORMAL( NCX , MCX , SCX ) 

KCY=NORMAL( NCY , MCY , SCY) 

KCP=NORMAI,(  NCP ,  MCP ,  SCP) 

*WIND  NOISE  (SLOW/FAST  VARYING) (VELOCITY  AND  DIRECTION) 

WVNS=NORMAL(  NWVS ,  UWVS ,  SWVS ) 

WVNF=NORMAL( NWVF , UWVF , SWVF ) 

WDNS=NORMAL( NWDS , UWDS , SWDS ) 

WDNF=NORMAL( NWDF , UWDF , SWDF) 

^ESTIMATOR  ERROR  -  INNOVATIONS  IN  SURGE,  SWAY,  AND  YAW 
INX=XT-XE 
INY=YT-YE 
INP=PT-PE 

’■'CURRENT  CORRECTION  (AS  APPLIED  TO  LF  VESSEL  MODEL) 

Q1=XD-UC 

Q2=YD-VC 

*THRUSTER  FORCE  =  BOW  THRUSTER  +  STERN  THRUSTER 
FT=BT+ST 

*BOW  THRUSTER  CONTROL 

BTI=  G21*EX+G22*EXXD+G23*£Y+G24*EYYD+G25*EP+G26*EPPD 
BTII=BTI-(DY*(VCH)*ABS(VCH)+FWYS) 

*  BTII=BTI-FWYS/2. 0 
BT=LIMIT( BTL, BTR , BTI I ) 

*STERN  THRUSTER  CONTROL 

STI=  G31*EX+G32*EXXD+G33*£Y+G34*EYYD+G35*EP+G36*EPPD 
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STII=STI-(DY*(VCH)*ABS(VCH)+FWYS) 

*  STII=STI-FWYS/2.0 
ST=LIM1T(STL,STR,STII) 

*PROP  CONTROL 

FPI=  GU*EX+G12*EXXD+G13*EY+G14*EYYD+G15*EP+G16*EPPD 
FPI I=FPI -( DX*UC*ABS( UCH)+0.  75*FWXS) 

*  FPII=FPI-FWXS 
FP=LIMIT(FPL,FPR,FPII) 

^MOMENT  CONTROL 

MT=78. 88*BT-81.  12*ST 

SAMPLE 

*TOTAL  VESSEL  MOTION  MODEL  -  MEASUREMENTS  IN  SURGE,  SWAY,  AND  YAW 
XT=X+XH+MNOIX 
YT=Y+YH+MNOIY 
PT*P+PH+MNOIP 

^SAMPLED  WIND  FORCES  AND  MOMENT  IN  SURGE,  SWAY,  AND  YAW 
FWXS=FWX 
FWYS=FWY 
MWS=MW 
DERIVATIVE 

.  *LOW  FREQUENCY  VESSEL  MODEL 
*LF  SURGE 

X=INTGRL(XO,XD) 

XD=INTGRL( XDO , XDD) 

XDD=NOIX+(FWX+FP-DX*Ql*ABS( Q1 )+MY*Q2*PD+MX*VC*PD) /MX 
*LF  SWAY 

Y=INTGRL(YO,YD) 

YD=INTGRL( YDO , YDD) 

YDD=NO I Y+( FWY+FT-DY*Q2*ABS ( Q2 )+MY*UC*PD+MX*Q 1*PD ) /MY 
*LF  YAW 

P=INTGRL( PO,PD) 

PD=INTGRL( PDO , PDD) 

PDD=NOIP+(-(MY-MX)*Ql*Q2-DP*ABS(PD)*PD+MC+MW+MT)/MP 
*1IIGH  FREQUENCY  VESSEL  MODEL 
’■'HF  SURGE 

XH=INTGRL(0. 0,XHD) 

XHD=INTGRL(0. 0,XHDD) 

XHDD=- ( WX**2 )*XH+HNOIX- . 2*XHD 
*HF  SWAY 

YH=INTGRL(0. 0,YHD) 

YHD=INTGRL(0. 0,YHDD) 

YHDD=-( WY**2 ) *YH+HNOI Y- . 2*YHD 
*HF  YAW 

PH=INTGRL( O.O, PHD) 

PHD=INTGRL( 0. 0 , PHDD) 

PHDD=- ( WP**2 )*PH+HNOIP- . 2* PHD 
*HF  WAVE  FREQUENCY  VARIATION  MODEL 
’•'FREQUENCY  IN  SURGE 

WX=INTGRL(W,WXD) 

WXD=NORMAL( NWX , UWX , SWX) 

’•'FREQUENCY  IN  SWAY 

WY=INTGRL(W,WYD) 

WYD=NORMAL( NWY , UWY , S W Y ) 

*FREQUENCY  IN  YAW 

WP=INTGRL(W,WPD) 


WPD=NORMAL( NWP , UWP , SWP) 

^CURRENT  MODEL  IN  EARTH  COORDINATES 
XCED=INTGRL(XCI ,XCEDD) 

XCEDD=CXNOI 

YCED=INTGRL(YCI ,YCEDD) 

YCEDD=CYNOI  . 

PCED=INTGRL( PCI , PCEDD) 

PCEDD=CPNOI 

’•'CURRENT  MODEL  IN  SURGE,  SWAY,  AND  YAW 
UC=XCED*DCOS( PT)+YCED*DSIN( PT) 
VC=-XCED*DSIN(PT)+YCED*DCOS(PT) 

MC=PCED 
*WIND  MODEL 

* SLOWLY  VARYING  WIND  VELOCITY 
WVS=INTGRL( IWVS , WVSD) 

WVSD=WVNS 

’•'RAPIDLY  VARYING  WIND  VELOCITY 
WVF=INTGRL(0. O.WVFD) 

WVFD=A1*WVF+WVNF 
*SLOWLY  VARYING  WIND  DIRECTION 
WDS=INTGRL( IWDS,WDSD) 

WDSD=WDNS 

’■'RAPIDLY  VARYING  WIND  DIRECTION 
WDF=INTGRL(0.  0,WDFD) 

WDFD=A2*WDF+WDNF 
’•'WIND  ANGLE  OF  ATTACK 
BETA=WDS+WDF-PT 

*WIND  FORCE  EQUATIONS  IN  SURGE,  SWAY,  AND  YAW 
FWX=-Fl’"'DCOS  ( BETA )  *(  WVS+WVF )  **2 
FWY=-F2*DSIN( BETA)*( WVS+WVF)**2 
MV=F3*DS IN( BETA )*( WVS+WVF) **2 
*KALHAN  FILTER  FOR  VESSEL 
’•'VESSEL  KALMAN  IN  SURGE 
EX=RX-XE 
EXXD=RXD-XED 

AXE=(-DX*(XED-0. )*ABS(XED-0.  )+MY*( YED-O.  )*PED+MX*0. *PED)/MX 

BXE=( FP+FWXS ) /MX 

XE=INTGRL(0.0,XED) 

XED=INTGRL(0. 0,XEDD)+K11*(XT-XE)+K12*(YT-YE)+K13*(PT-PE) 
XEDDaAyS+BXE+K2 1*( XT-XE)+K22*( YT-YE )+K23*( PT-PE ) 

* VESSEL  KALMAN  IN  SWAY 
EY=RY-YE 
EYYD=RYD-YED 

AYE=( -DY*( YED-O. )*ABS(YED-0.  )+MX*(XED-0.  )*PED+MY*0.  *PED)/MY 

BYE=(FT+FWYS)/MY 

YE=INTGRL(0.  0,YED) 

YED=INTGRL(0. 0,YEDD)+K31*(XT-XE)+K32*(YT-YE)+K33*(PT-PE) 
YEDD=AYE+BYE+K41*(XT-XE)+K42*(YT-YE)+K43*(PT-PE) 

^VESSEL  KALMAN  IN  YAW 
EP=RP-PE 
EPPD=RPD-PED 

APE=(-(MY-MX)*(XED-0.  )*(YED-0.  )-DP*PED*ABS(PED))/MP 
BPE=(MWS+MT+0. )/MP 
PE=INTGRL( O.O, PED) 

PED=INTGRL(0. 0,PEDD)+K51*(XT-XE)+K52*(YT-YE)+K53*(PT-PE) 


PEDD=APE+BPE+K61*(XT-XE)+K62*(YT-YE)+K63*(PT-PE) 

*FILTER  PREDICTION  OF  CURRENT  OFFSET 
^CURRENT  PREDICTION  IN  SURGE 
CX=INX-XCH 
XCH=INTGRL(0.  0,UCH) 

UCH=INTGRL(0.  0,UCHD)+(C11*CX) 

UCHD=(C21*INX) 

^CURRENT  PREDICTION  IN  SWAY 
CY=INY*YCH 
YCH=INTGRL(0.  0,VCH) 

VCH=INTGRL(0.  0,VCHD)+(C32*CY) 

VCHD=(C42*INY) 

’•'CURRENT  PREDICTION  IN  YAW 
CP=INP-PCH 
PCH=INTGRL(0.  0,PCHD) 

PCHD=INTGRL( 0. 0 , PCHDD)+( C53*CP) 

PCHDD=(C63*CP) 

*  TRANSFORM  BACK  TO  VESSEL  COORDINATES 
MCH-PCHD 

CONTROL  FINTIM=800. 0,DELT=. 10,DELS=1.0 
METHOD  RKSFX 

SAVE  10.  OOjX^XjFPjYjEYjBTjSTjMT^jEP.XTjYTjPTjXEjYEjPEjXHjYI^PH,.  . 

FWX , FWY , MW , XDD , XHDD , MC , UC , VC , UCH , VCH , INX , INY , INP , XED , YED , PED 
*RINT  UC, UCH, VC, VCH 

GRAPH  ( G 1 ,DE=TEK618)  TIME(UN=' SEC' ) ,XT(UN=' FT’ ) ,XE 
LABEL  (Gl)  XT  VS  TIME 

GRAPH  (G2,DE=TEK618)  TIME(UN=' SEC' ) , YT(UN=' FT' ) ,YE 
LABEL  (G2)  YT  VS  TIME 

GRAPH  (G3,DE=TEK618)  TIME(UN=' SEC* ) ,PT(UN='RAD' ) ,PE 
LABEL  (G3)  PT  VS  TIME 

GRAPH  (G4,DE=TEK618)  TIME(UN='SEC' ) ,FP(UN='LBF' ) 

LABEL  (G4)  PROP  CONTROL 

GRAPH  (G5,DE=TEK618)  TIME, BT, ST 

LABEL  (G5)  THRUSTER  CONTROL 

GRAPH  (G6,DE=TEK618)  Y(UN='FT' ) ,X(UN='FT' ) 

LABEL  (G6)  POSITION  X  VS  Y 
GRAPH  (G7,DE=TEK618)  XE,XED 
LABEL  (G7 )  PHASE  PLANE 
*GRAPH  (G7,DE=TEK618)  TIME, UC, UCH 
*LABEL  (G7)  POSITION  UC,UCH  VS  TIME 
*GRAPH  (G8,DE=TEK618)  TIME, VC, VCH 
*LABEL  (G8)  POSITION  VC, VCH  VS  TIME 
*GRAPH  (G9,DE=TEK618)  TIME, INX, INY 
’’'LABEL  (G9)  ERROR  RESIDUALS 
’•'GRAPH  (G8,DE=TEK618)  TIME,Y,YH 
’•'LABEL  (G8)  HIGH  AND  LOW  FREQUENCY  SWAY 
’•'GRAPH  (G9,DE=TEK618)  TIME,P,PH 
’•'LABEL  (G9)  HIGH  AND  LOW  FREQUENCY  YAW 
END 
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